Preprint typeset in JHBP style - HYPER VERSION 



OU-HEP-111006 



Coupled Boltzmann calculation of 
mixed axion/neutralino cold dark matter 
production in the early universe 



Howard Baer", Andre Lessa'' and Warintorn Sreethawong" 

'^Dept. of Physics and Astronomy, University of Oklahoma, Norman, OK 73019, USA 
^ Instituto de Fisica, Universidade de Sao Paulo, Sao Paulo - SP, Brazil 
E-mail: 



baerOnhn . ou . edu 



lessaOf ma. if .usp.br, wstanOrLliii. ou. edu 



Abstract: We calculate the relic abundance of mixed axion/neutralino cold dark matter 
which arises in i2-parity conserving supersymmetric (SUSY) models wherein the strong 
CP problem is solved by the Peccei-Quinn (PQ) mechanism with a concommitant ax- 
ion/saxion/axino supermultiplet. By numerically solving the coupled Boltzmann equa- 
tions, we include the combined effects of 1. thermal axino production with cascade decays 
to a neutralino LSP, 2. thermal saxion production and production via coherent oscillations 
along with cascade decays and entropy injection, 3. thermal neutralino production and 
re-annihilation after both axino and saxion decays, 4. gravitino production and decay and 
5. axion production both thermally and via oscillations. For SUSY models with too high a 
standard neutralino thermal abundance, we find the combined effect of SUSY PQ particles 
is not enough to lower the neutralino abundance down to its measured value, while at 
the same time respecting bounds on late-decaying neutral particles from BBN. However, 
models with a standard neutralino underabundance can now be allowed with either neu- 
tralino or axion domination of dark matter, and furthermore, these models can allow the 
PQ breaking scale fa to be pushed up into the 10^^ — 10^^ GeV range, which is where it is 
typically expected to be in string theory models. 
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1. Introduction 



The Standard Model (SM) of particle physics is beset by two afflictions: 1. in the scalar 
(Higgs) sector of the theory, quadratic divergences require large fine-tunings of electroweak 
parameters which depend on the scale A below which the SM is regarded as the correct 
effective field theory of nature and 2. in the QCD sector of the theory, the Lagrangian 
term 

^ 5 4^2PA,.F'i^ (1-1) 

required by 'tHooft's solution to the U{1)a problem is constrained to a value 9 < 10^^*^ 
to gain accord with measurements of the neutron EDMjl|. The first of these is solved 
by the introduction of softly broken weak scale supersymmetry (SUSY) into the theory||2[ 
(which receives some indirect support from the measured values of gauge couplings at 
LEP|^] and from global fits to precision electroweak data|Q]), while the second problem 
is solved by the introduction of a global U{1)pq Peccei-Quinn (PQ) symmetry broken by 
QCD anomalies which requires the existence of an ("invisible") axion|^, ^, with mass 
expected in the micro-eV or below range ||8|. Solving both problems simultaneously requires 
supersymmetrization of the SM (usually via the Minimal Supersymmetric Standard Model, 
or MSSM) along with the introduction of an axion supermultiplet a into the theory. The a 
supermultiplet contains an i?-parity-even spin-0 saxion field s{x) along with an i?-parity- 
odd spin-^ axino a{x), in addition to the usual pseudoscalar axion field a{x): 

a = + iV26aL + iOeL^a, (1.2) 

v2 

in 4-component spinor notation[^. 

In such a theory, it is expected that SM superpartner particles with weak scale masses 
should emerge, along with a weak scale saxion, whilst the axino mass is more model de- 
pendent, with rria ~ keV-TeV being expected]^]. The axion, saxion and axino couplings 
to matter depend on the PQ breaking scale fa^, which is required fa > 10^ GeV by stellar 
cooling calculations |llO|] . The axion is often considered as a very appealing dark matter 



(DM) candidate 111, 



In the MSSM, DM candidates include the lightest neutralino Zi (a WIMP), the spin-| 
gravitino G or possibly the superpartner of a right-handed neutrino |T^ . Gravitino dark 
matter is tightly constrained and disfavored by the standard picture of Big Bang nucle- 



osynthesis (BBN)|15], whilst right-hand neutrino states are expected to exist near the GUT 
scale according to the elegant see-saw mechanism for neutrino mass 10]. Many authors thus 
expect dark matter to be comprised of the SUSY neutralinos, a natural WIMP candidate 
which is motivated by the so-called "WIMP miracle". However, detailed analyses show 
that neutralino dark matter requires a rather high degree of fine-tuning[17| to match the 



^Throughout this work we omit the number of generations factor A'^, which appear along with the PQ 
scale, fa/N, in the DSFZ model and in the KSVZ model with more than one heavy quark generation. All 
our results can then be trivially generalized replacing fa by fa/N. 

^For a somewhat different axion/axino scenario, see Ref. Iisl. 
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WMAP-measured cold DM abundance fl^: 

Qj^uh^ = 0.1123 ± 0.0035 at 68% CL. (1.3) 

In fact, the measured abundance lies in the most improbable locus of values of neutralino 
relic density as predicted by general scans over SUSY model parameter space []l9| . 

The PQ-extended Minimal Supersymmetric Standard Model (PQMSSM) offers addi- 
tional possibilities to describe the dark matter content of the universe. In the PQMSSM, 



the axino may play the role of stable lightest SUSY partner (LSP)pO|, 21 1, while the quasi- 
stable axion may also constitute a component of DM|2^, giving rise to mixed axion/axino 
{ad) CDM. In supergravity theories however, the axino mass is expected to lie at the 
weak scale [p^], so that the neutralino remains as LSP, and the possibility occurs for mixed 
axion/neutralino (aZi) CDM. 



In a recent paper, Choi et al.\24] presented a semi-analytic approach for estimating 
the relic abundance of neutralinos in the mixed aZi CDM scenario. This approach applies 
to cases where the thermally averaged neutralino annihilation cross section times relative 
velocity (av) is approximately constant with temperature, as occurs for a wino-like or 



higgsino-like neutralino [ 25 1. Detailed calculations of the relic abundance of mixed aZi 
CDM were performed in Ref. |26|, where formulae for the neutralino and axion abundances 
were presented. 

The standard calculation of the neutralino Yield Y-''^ = (where n^^ is the neu- 
tralino number density and s is the entropy density) gives 

Zi 4{av)MpTfr ' ^ ' ' 

where g*{Tfr) is the number of active degrees of freedom at temperature T = Tfr, where 

3^/5(crw)Mpm~/^ 
Tf/ = m,Jln[ ^M . (1.5) 

T^^/^Tj^. g* [Tfr) 

is the freeze-out temperature and Mp is the reduced Planck mass. 

If instead axinos are thermally produced (TP) at a large rate at re-heat temperature 
Tr after inflation, then they cascade decay to (stable) neutralinos at decay temperature 

n = v/rVWp/ (7r25,(rB)/90)'/', (1.6) 

and can boost the neutralino abundance. The late-time injection of neutralinos into the 
cosmic soup at temperatures < Tfr may cause a neutralino re- annihilation effect such 
that the neutralino Yield is instead given by|2^, ^ 

Zi '^-^S 4{av)MpT^ ■ ^ ' ^ 

Since is typically in the MeV-GeV range, i.e. well below Tfr ~ the neutralino 

abundance after re-annihilation can be highly enhanced relative to the standard cosmo- 
logical picture. In addition, one must fold into the relic abundance the axion contribution 
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arising from coherent axion field oscillations beginning at axion oscillation temperature 
Ta ~ 1 GeV. 

An additional complication comes from entropy production from axino decay after 
Tfj. (which may dilute the neutralino abundance) or after Ta (which may dilute the axion 
abundance). This may occur in the case where axinos temporarily dominate the energy 
density of the universe. Depending on the PQ parameters of the PQMSSM model (/«, m^, 
initial axion misalignment angle Oi and Tr), the dark matter abundance may be either 
neutralino- or axion-dominated. In fact, cases may occur where the DM relic abundance 
is shared comparably between the two. In the latter case, it might be possible to directly 
detect relic neutralino WIMP particles as well as relic axions! 



While the semi-analytic treatment of Ref's |24| and [^] provides a broad portrait of 
the mixed aZ\ CDM picture, a number of important features have been neglected. These 
include the following. 

• For bino-like neutralinos, {av) ~ a -|- where a ~ since we mainly have p-wave 
annihilation cross sections. In this case, {av) is no longer independent of temperature. 



and the simple formulae 1.4 and 1.7 are no longer valid. 



• In Ref's ||2j] and [g6|], the effects of saxion production and decay in the early universe 
are neglected. In fact, saxion thermal production or production via coherent oscilla- 
tions (C0)[|^], followed by late time saxion decay, may inject considerable entropy 
into the early universe, thus diluting all relics present at the saxion decay temper- 
ature T^. Saxions may also add to the neutralino abundance via decays such as 
s — )• gg, followed by gluino cascade decays. There exists the possibility of saxion and 
axino co- domination of the universe. In this case, there might be a second neutralino 
re-annihilation taking place at Tf-^. 



The treatments of ||2J] and |2^] invoke the "sudden decay" approximation for late- 
decaying axinos, whereas in fact the decay process is a continuous one proceeding in 
time until the decaying species is highly depleted (all have decayed). 



• The treatments of [^4[ and [26| largely ignore the effect of gravitino production and 
decay in the early universe. 

To include the above effects into a calculation of the mixed aZi relic abundance, one 



must go beyond the semi-analytic treatment presented in Ref's |24, and proceed with 
a full solution of the coupled Boltzmann equations which govern various abundances of 
neutralinos, axinos, axions, saxions, gravitinos and radiation. 

Toward this end, in Sec. ^ we present a simplified set of coupled Boltzmann equations, 
which we use to calculate the relic abundance of mixed axion/neutralino dark matter. 
More details about the approximations made and each term present in our equations are 
discussed in Appendix ^ 

In Sec. 1^, we present various numerical results for the mixed aZi CDM scenario 
using the full set of Boltzmann equations. We find that, even after the inclusion of the 
saxion field, adjusting the parameters of the PQMSSM can only increase the neutralino 
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abundance, and not decrease it, while at the same time respecting bounds on late-decaying 



neutral particles from BBN. This result is the same as found in Refs. |24| and |26|, but 
now corresponds to a much stronger statement, since the saxion entropy injection had 
been neglected in the previous works. Furthermore, our results also apply to models with 
bino-like neutralinos, which could not be studied in the semi- analytical framework used in 
Refs. H and [M. 



Since the neutralino abundance can be only enhanced in the PQMSSM, in models such 
as mSUGRA, those points which are excluded by a standard overabundance of neutralinos 
are still excluded in the PQMSSM! This rather strong conclusion does depend on at least 
three assumptions: 1. that thermal axino production rates are not suppressed by low-lying 



PQ-charged matter multiplets[28p, 2. that saxion decay is dominated by gluon and gluino 
pairs and 3. that the assumed saxion field strength s(x) = 9sfa is of order the PQ-breaking 
scale fa, i.e. that 0^ ~ 1- 

We also examine several cases with a standard underabundance of neutralino dark 
matter. In these cases, again the neutralino abundance is only increased (if BBN constraints 
are respected). Thus, adjustment of PQMSSM parameters can bring models with an 
underabundance of neutralinos into accord with the measured DM relic density. In these 
cases, the DM abundance tends to be neutralino-dominated. Also, in these cases, solutions 
exist where the PQ scale fa is either near its lower range, or where fa is much closer to 
Mgut, with fa ~ 10^^ GeV typically allowed. This is much closer to the scale of fa which is 
thought to arise from string theory|31|. In Sec. ^, we present a summary and conclusions. 



2. Mixed axion/neutralino abundance from coupled Boltzmann equations 

Here, we present a brief description of our procedure to calculate the relic abundance of 
mixed aZi CDM in the PQMSSM. A more detailed discussion is left to Appendix 

2.1 Boltzmann equations 

The general Boltzmann equation for the number density of a particle species can be gener- 



ically written as|32|: 



fii + ^Hui = Si - —Tim (2.1) 

It 

where Si represents a source term, Fj is the decay width and 7^ is the relativistic dilation 
factor to take into account the suppressed decays of relativistic particles. To describe the 
thermal production of a particle species i as well as its decoupling from the radiation fluid 
and the non-thermal production coming from other particles decays, we include in Si the 
following terms: 

S. = -\n\ - {nf{T))\avUT)^Y.^m^)^o- (2-2) 



^Here, we assume standard rates for thermal axino production as calculated in the literature j2^, po[ . 
In Ref. [Q, it has been shown that if PQ-charged matter multiplets $ exist well below the PQ breaking 
scale /a, then axino production is suppressed by factors of vn^jTR. 
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where {av) is the (temperature dependent) thermally averaged annihilation cross section 
times velocity for the particle species z, n^"^ is its equilibrium number density and BR{j, i) 
is the branching fraction for particle j to decay to particle i.^ 
The Boltzmann equation then becomes: 



n. 



+ 3Hn, = -Fmi— + [«'(r))2 - n^,]{av)i + ^ 5i?(j, i)r,m,-^ , (2.3) 

Pi j Pj 



where we have used 7^ = pi/miTii. As discussed in Appendix ^ the above equation is also 
valid for coherent oscillating fields once we take BR{j, i) = and {av)i = 0. 
It is also convenient to write an equation for the evolution of entropy: 

.2 1 X 1/3 



or S = —Y,BRii,X)rimm (2.4) 

i 

where BR{i,X) is the fraction of energy injected in the thermal bath from i decays. 
Along with Friedmann's equation, 



the set of coupled differential equations, Eq's. |2.3| , p.4| and |2.5| , can be solved as a function 
of time. More details on the solution of the above equations and the expressions used for 
{av)i, BR{i,j) and BR{i,X) are presented in Appendix^. 

2.2 Present day abundances and constraints from BBN 

To compute the relic density of neutralinos and axions we evolve the various particle and 
sparticle abundances from T = until the final temperature Tp is reached at which all 
unstable particles (save the axion itself) have decayed. The relic densities of the various 
dark matter species labeled by i are then given by: 

^.h^ = ^ X (2.6) 
s{Tf) Pc/h^ 

In our calculations, a critical constraint comes from maintaining the success of the 
standard picture of Big Bang nucleosynthesis. Constraints from BBN on late decaying 



neutral particles (labeled X) have been calculated recently by several groups fSj, |36| 



(we explicitly use the results of Ref. [^]) and are presented as functions of 1. the decaying 



*In this paper, i is summed over 1. neutralinos Zi, 2. TP axinos a, 3. and 4. CO- and TP-produced 
saxions s{x), 5. and 6. CO- and TP-axions a, 7. TP gravitinos G and radiation. We allow for axino 
decay to gg, ^Zi and ZZi states {i — 1 — 4), and saxion decay to gg and gg. Additional model-dependent 
saxion decays e.g. to aa and/or hh are possible and would modify our results. We assume G decay to all 
particle-sparticle pairs, and include 3-body gravitino modes as well|33]. 
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neutral particle's hadronic branching fraction B^, 2. the decaying particle's lifetime tx, 
and 3. the decaying particle's relic abundance Qxh'^ had it not decayed. The constraints 
also depend on 4. the decaying particle's mass mx- We have constructed digitized fits 
to the constraints given in Ref. |^^, and apply these to late decaying gravitinos, axinos 
and saxions. Typically, unstable neutrals with decay temperature below 5 MeV (decaying 
during or after BBN) and/or large abundances will be more likely to destroy the predicted 
light element abundances. 

2.3 Example: calculation from a generic mSUGRA point 

As an example calculation, we adopt a benchmark point from the paradigm minimal super- 
gravity model (mSUGRA), with parameters (mo, mi/2, Aq, tan/3, sign{fi)) = (400 GeV, 
400 GeV, 0, 10, +). The sparticle mass spectrum is generated by Isasugrap^, and has 



a bino-like neutralino with mass uit^ = 162.9 GeV and a standard relic abundance from 
IsaReD|38] of Q.^-'^h? = 1.9 (it would thus be excluded by WMAP7 measurements assuming 



the standard neutralino freeze-out calculation). We assume a gravitino mass = 1 TeV. 

Here, we work in the PQMSSM framework, and take Tr = 10^^ GeV with PQ param- 
eters as ma = 1 TeV, m^ = 5 TeV, Oi = 0.5 and fa = 10^^ QgV. We also take Os = 1, 
where 9s fa is the initial field amplitude for coherent oscillating saxions. The various energy 
densities pi are shown in Fig. || for i = R (radiation), Zi (neutralinos), a^^ (thermally 
produced axions), a^^ (coherent oscillating axions), s^^ (thermally produced saxions), 
s^'-^ (coherent oscillating saxions), a^^ (thermally produced axinos) and G^^ (thermally 
produced gravitinos). The energy densities are plotted against scale factor ratio R/Rq, 
where Rq is the scale factor at T = Tr. We also plot the temperature T of radiation 
(green-dashed curve). 

We see that, at R/Rq < 10^°, the universe is indeed radiation-dominated. At T ^ 1 
TeV, the TP axions, saxions and axinos all have similar abundances. At these temperatures, 
the saxion coherent abundance as well as the gravitino thermal abundance are far below the 
other components. As the universe expands and cools, most components are relativistic, 
and decrease with the same slope as radiation: pi ~ T"^. The exception is the CO- 
produced saxions, which are non-relativistic, and fall-off as p^*-^ ~ T~^. At R/Rq ~ 10'', the 
temperature T ~ 1 TeV, and the thermally-produced axinos, saxions and gravitinos become 
non-relativistic, so now p^^ ~ ~ For even lower temperatures with R/Rq ~ 10 , 

neutralinos begin to freeze-out, and their abundance falls steeply. At T ~ m^^/20, they do 
freeze-out, and normally their density would fall as T^^, as indicated by the blue dot- 

dashed curve, which shows neutralino abundance in the MSSM, without PQ- augmentation. 
In the PQMSSM however, saxions- and later still axinos- begin decaying in earnest, and 
feed into the neutralino abundance, preventing its usual fall as T~^. At T ~ 0.5 GeV, 
the energy density of axinos surpass the radiation component and the universe becomes 
axino-dominated until the axino decays at T ~ 10 MeV. Also, around R/Rq ~ 3x 10^° with 
T ~ 1 GeV, CO production of axions begins, and by R/Rq ~ 4 x 10^^, with T < Aqcd, its 
abundance begins to fall as T-3. For even lower temperatures (T < 10 MeV), the axinos 
have essentially all decayed, feeding back into the neutralino abundance, and also increasing 
the entropy per co- moving volume, which would otherwise be conserved. At R/Rq ~ 10^^, 
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Figure 1: Evolution of radiation, neutralino, axion, saxion, axino and gravitino energy 
densities versus scale factor R. We adopt an mSUGRA SUSY model with parameters 
(mo,mi/2, Ao,tan^, si5n(/x)) = (400 GeV,400 GeV,0, 10,+). We also take = 1 TeV and 
Tr = 10^° GeV and PQ parameters m~a = 1 TeV, = 5 TeV, = 0.5, 9^ = 1 with fa = lO^^ 
GeV. 



the universe moves from radiation domination to matter (neutralino) domination, while 
at even lower temperatures, the gravitinos decay away. In this case, the final neutralino 
abundance is ^^^^^^ ~ 90017- far beyond its standard value. This is mainly due to its 
abundance being augmented by thermal axino and saxion production and cascade decay 
to neutralinos. In the standard axion cosmology, the axion abundance would have been 
Qstdj^2 ^ 0.06p9[|. In the case illustrated here, entropy injection from saxion, axino and 
gravitino decays has diluted its abundance to just fia^^ ~ 0.004. 

As an example of the relevance of using the full set of Boltzmann equations instead of 



the semi- analytical approach of Refs. [24| and [^^, we compare in Fig. g the neutralino and 
axion relic densities as a function of the axino mass using the Boltzmann equation formalism 
and the semi-analytical approach. The other PQMSSM parameters are the same as used 
in Fig. H], but to compare with the semi-analytical results of Refs. [p4| and |26] we neglect 
the saxion component. For these choices of PQ parameters and for rria ^ 50 TeV, the axino 
decays after neutralino freeze-out (as seen on Fig. |l|) , significantly enhancing its final relic 
abundance. Furthermore, the axino decay injects entropy, diluting the axion abundance. 
As we can see from Fig. the axion relic density obtained using the analytical expressions 
derived in Ref. [p6| agree extremely well with the solution of the Boltzmann equations. On 
the other hand, the analytic neutralino abundance disagrees with the Boltzmann solution 
by almost an order of magnitude for rria ^ 50 TeV. The primary reason for this is the fact 
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Figure 2: Neutralino and axion relic densities as a function of the axino mass for Qi — 0.5, 
Tr = 10^° GeV, fa = 10^2 GeV and the mSUGRA point (toq, mi/2, Aq, tan/3, si5n(/i)) = 
(400 GeV, 400 GeV, 0, 10, +). The sohd hnes correspond to the solution of the Bohzmann equations 
while the dashed lines correspond to the results obtained using the analytical expressions derived 
inRef. [§6l. 



that, for this mSUGRA point, the neutralino is bino-like and {av)^^ is no longer constant, 
but strongly depends on the temperature. This dependence has not been included in the 
semi-analytical approach. Also, the sharp transition seen in the semi-analytical result 
at rria — 18 TeV, where becomes bigger than Tf^, is artificially introduced by the 
sudden decay approximation. As shown by the Boltzmann solution, the enhancement of 
the neutralino relic abundance smoothly decreases, going up to ~ 50 TeV. 

3. Neutralino abundance in the PQMSSM 

3.1 Neutralino Abundance in several PQMSSM models 

In this section, we adopt four SUSY benchmark models listed in Table ||. The first two 
points, labeled BMl and BM2, are generic mSUGRA points with a bino-like Zi which 
give rise as expected to an apparent excess of CDM. For BMl, with (mo, mi/2j ^O; tan/3, 
sign{n)) = (400 GeV, 400 GeV, 0, 10, +) we have m^^ = 162.8 GeV with a standard abun- 
dance il'-'^h? = 1.9, while the second point (BM2) has (mo, ^-1/2; ^O; tan/3, sign{^)) = 
(3000 GeV, 1000 GeV, 0, 10, +) with m^^ = 436.3 GeV and Vt'i'^^h^ = 49.6. The next 
point, BM3, is a mSUGRA point with an apparent underabundance of neutralino dark 
matter, with = 163.8 GeV, tua = 367.5 GeV, lying in the ^-funnel region|Q, so 
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BMl 


BM2 


BM3 


BM4 


mo 


400 


3000 


400 









1 nnn 






"13/2 


10^ 


10^ 


10^ 


5 X 10^ 


tan /3 


10 


10 


55 


10 




162.9 


436.3 


163.8 


142.1 


a^^{Zip) pb 


1.9 


49.6 


0.019 


0.0016 


8.1 X 10-1° 


1.1 X 10-1° 


2.1 X 10-s 


4.3 X 10-9 



Table 1: Masses and parameters in GeV units for several benchmark points computed with 
Isajet7.81 using = and nit = 173.3 GeV. 



^stdf^2 _ Q QiQ Pqj. a,ll cases, we take = 1 TeV, but now will vary the PQ parameters 
and Tr, in order to see if the relic density of mixed aZi CDM can lie in the WMAP-allowed 



region. The last point is taken from the gaugino AMSB model[41, 42] and has a wino-like 
neutralino with Q^^'^h'^ = 0.0016, but with = mg = 50 TeV. 

In order to keep our results as general as possible, we will not assume particular PQ 
parameters, but instead we scan over the following parameter values: 

10^ GeV < fa < 10^^ GeV, (3.1) 

500 GeV <ma < 10^ GeV, (3.2) 

10^ GeV < ms < 10^ GeV, (3.3) 

0.1 < Os < 10, (3.4) 

10^ GeV <Tr < 10^2 GeV. (3.5) 

Since we will be mostly concerned with the neutralino relic abundance, we leave the axion 
mis-alignment angle 9i undetermined for now. 

3.2 Benchmark BMl 

Our results are shown as the resultant relic density of neutralinos ^^^h'^ in the PQMSSM, 
where we plot each model versus fa in Fig. ^. The blue points are labeled as BBN-allowed, 
while red points violate the BBN bounds as described in Sec. |2.2| . 

From Fig. ^, we see that at low values of PQ breaking scale fa ~ 10^ — lO^i GeV, the 
value of ^^_^h'^ is always bounded from below by its standard value Q^-^h? ~ 1.9. Those 
points with ^^^h? ~ 1.9 are typically those for which axinos and saxions decay before Tf^,, 
or those for which axino/saxion production is suppressed by low Tr so that axinos/saxions 
decays do not significantly contribute to Q.^Ji'^. 

In Fig. 1^ frequently the neutralino abundance is enhanced beyond 1.9, making these 
points even more excluded. The reason why points only have enhanced relic densities at 
the lower fa range is because the axino-matter coupling is proportional to l//a, and so 
thermal axino production is enhanced compared to higher fa values. In addition, a ^ gg 
decays may be phase space suppressed, so that axino decay takes place at temperature 
Tf) < Tfr, thereby augmenting the neutralino abundance. Saxion decay is never phase 
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Figure 3: Calculated neutralino relic abundance from mSUGRA model BMl versus fa- We take 
itIq = 1 TeV. The spread in dots is due to a scan over PQ parameters fa, Tji, rria, tris and Og- 



space suppressed, since s ^ gg is always possible, so at the lower range of fa, saxion decay 
typically takes place at > Tjr. For values of fa ~ 10^^ GeV and beyond, axinos can 
no longer decay before neutralino freeze-out, and so the neutralino abundance is always 
enhanced. The value of fa where the neutralino abundance is always enhanced is somewhat 
an artifact of our scanning range, since if we allow > 10^ GeV, axinos could become 
shorter- lived for a higher value of fa since Fq ~ fa- 

At even higher values of fa ^ 10^^ GeV, axino/saxion thermal production becomes 
increasingly suppressed, while saxion production via CO becomes enhanced: entropy dilu- 
tion by saxions starts winning over neutralino production from thermal axino production 
and decay. Also, both saxion and axino become even longer-lived, so more points become 
BBN-disallowed. Although the entropy injection from s ^ gg decays grows with fa, the 
BBN-allowed blue points are never pushed below ^^x^"^ ^ s™ce s gg also injects 
additional neutralinos into the thermal bath.^ 

To understand why the neutralino injection from saxion decays always wins over the 
entropy dilution of the neutralino abundance, we must look at the neutralino Yield from 
saxion decays. For simplicity, we will neglect the axino component as well as neutralino 
re-annihilation at Tfy and assume that the PQ parameters are chosen so the universe has 
a saxion-dominated era, since this is the only scenario with significant entropy injection. 

^We checked the effect of artificially turning off s — s- decays in Fig. ^. In this case, at high fa > 10^'' 
GeV, the enhanced saxion production via COs produces only entropy dilution of the neutralino abundance, 
and some BBN-allowed points remain with a highly suppressed neutralino abundance at high fa . This effect 
has lead to claims that large fa ~ Mqut values may be allowed in SUSY models due to entropy injection 
by saxions^, ^ |50[|l|, |52j . By properly including s gg decay and the concommitant 

neutralino re-annihilation at Ts, the pure entropy injection effect is counter-balanced in this case, and the 
highly diluted cases become BBN-forbidden. 
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Under these assumptions, the Yield of neutralinos emitted from saxion decays is simply 
given by 1 52]: 

Y^^ = h^s X 2BR{s ^ g-g) (3.6) 

where the factor 2 above takes care of the multiplicity of neutralinos from each saxion 
decay and r is the entropy injection factor, which can be approximated by 

r = ^. (3.7) 



where Tf^ is the saxion decay temperature and Tg = AmgYg/'i (see Refs. [g2|, |5§]). Therefore: 



= -^BR(s gg) =^ n% 4 x 10^ GeV"! '^BR(s gg). (3.8) 
'^i 2 iris ^ 'Tig 

The above expression shows that the relic density can be suppressed for large nis , small Tf-^ 
and/or small BR. However, as seen in Fig. |3|, such suppression never seems to drive ^^^h"^ 
below its standard value, except in the BBN excluded region. To see why this happens, 
using Eq. 3.8 we show in Fig. contours of Q% h? in the vs. fa plane for the BMl point. 



with Tr = 10*^ GeV and Os = I. We also show the BBN excluded region (Tf, < 5 MeV) 
and the region with Tf^ > Tg (r < 1), where there is no saxion dominated era and Eq. |3.8| 
is no longer valid. As we can see, the allowed region (white) can only satisfy the WMAP 
constraints at very large iJis and /„ values. The main reason for the low ^^J^^ values 
obtained in this region is due to the suppression of BR{s — )• gg). This can be seen in Fig. 
|5|, where we show the branching ratio as a function of nig for the same benchmark point. 
We can see that- for the region where the s ^ gg decay mode is closed- the saxion lifetime 
falls into the BBN-forbidden zone. This can also be seen in Fig. ^, where all the low ^^_^h? 
points at large fa are BBN excluded. 

As seen from the above results, the neutralino relic abundance can indeed be diluted 
by including the saxion field, but only at the expense of going to extremely high and 
fa values. However, since is expected to be of order the soft SUSY masses (or ~ in 
e.g. AMSB models), we consider such high values extremely unnatural. Furthermore, the 
PQMSSM is most likely not the correct effective theory at fa > 10^^ — 10^^ GeV, where 
we expect a Grand Unified theory and/or large supergravity corrections. Nonetheless, to 
confirm the approximate results obtained from Eq. |3.8| , we extend our previous scan over 
to 

fa G [10l^ 10^2] GeV , rus G [10^ 10^] GeV (3.9) 

and use the full set of Boltzmann equations to compute the neutralino relic abundance. 
The results are shown in Fig. ^, where we plot in the vs. fa plane all solutions satisfying 
l^^^/i^ < 0.11. As we can see, the numerical results agree very well with the analytical 
results in Fig. ^. The only discrepancy is in the region near = T^, which does not 
present viable solutions in the scan. This is simply due to the fact that in our estimate of 



Eq. 3^, we neglected the neutralino freeze-out component, which becomes dominant when 
r ~ 1 or Tg ~ T^, increasing the value of ^^_^h'^ in this region. 
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m~ = 163 GeV, m = 951 GeV, r„ = 10^ GeV, 6=7 




10^ 10* 10^ 10'' 10^ 

rris (GeV) 

Figure 4: Regions in the iris vs. fa plane where T|) < 5 MeV and T|) > Tg. We also show contours 



of constant fi- h as estimated using Eq. 3.8 

Zi I 1 



From the results presented above, we see that for reasonable values of fa, rUg and Og, 
the result illustrated for point BMl in Fig. ^ seems to generalize to all SUSY model points 
with a standard neutralino overabundance: SUSY models with a standard overabundance of 
neutralino dark matter are still at least as excluded when augmented by the PQ mechanism. 

In Fig. ^, we show the evolution of various energy densities versus the scale factor 
for a large fa value. In this case, we see that the universe is radiation-dominated out 
to R/Ro ~ 10^, whereupon it becomes saxion dominated. If only s ^ gg is considered, 
the saxion entropy injection would cause a large dilution of neutralinos. But by including 
s ^ gg decays, we see the neutralino enhancement during 10'' < R/Rq < 10^. We also 
show by the dash-dotted line the expected neutralino energy density by neglecting s ^ gg 
decays: in this case, the neutralino abundance is highly suppressed compared to the case 
where s ^ gg \s accounted for. 

3.3 Benchmark BM2 

To emphasize some of the generality of our previous results, we show a further point with 
a standard overabundance of neutralinos in Fig. ^, with (mo, ^O; tan /3, si(7n(^)) = 
(3000 GeV, 1000 GeV,0, 10,+), for which VL%'^h^ ~ 50. By scanning over PQ parameters, 
again we find that for low /a, '^^^^ either remains at its standard value (if axinos/saxions 
decay before freeze-out), or are enhanced (if axinos/saxions decay after freeze-out). At 
high /a, entropy dilution from CO-produced saxions again can suppress /i^, but the 
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suppression is counterbalanced hy s ^ gg decays: only BBN-forbidden points where nis is 
so light that s ^ gg is kinematically closed yield points with ^^_^h'^ < 0.11. 

3.4 Benchmark BM3: A-resonance region 

In Fig. ^, we show the neutralino abundance in the case of an mSUGRA point lying in the 
yl-resonance annihilation region [pO|] where 2m^ ~ ttia- We adopt mSUGRA parameters 
{mo, mi/2,Ao, tan p,sign{fi)) = (400 GeV,400 GeV,0,55,+), for which Q'^'^h^ ~ 0.02, i.e. 
a standard underabundance.^ In this case, a scan over PQ parameters yields many points 
at low fa with ^^^h? ~ 0.02 — 10. Thus, the standard neutralino underabundance may 
be enhanced up to the WMAP-allowed value, or even beyond. As we push to higher fa 
values, the axino becomes so long-lived that it only decays after neutralino freeze-out, and 
hence the neutralino abundance is always enhanced. Above fa ~ 10^^ GeV, the neutralino 
abundance is enhanced into the WMAP-forbidden region, with ^^^Ji^ always larger than 
0.11. As we push even higher in /„, then axino production is suppressed, but CO-production 
of saxions becomes large. Entropy dilution turns the range of ^^^^^^ back down again, and 
at fa ~ 10^^ GeV, some BBN-allowed points again reach ^^^h? ~ 0.11. In this case, rather 
large fa values approaching Mqut are allowed. 

For the points with ^^Ji'^ < 0.11, the remaining dark matter abundance can be 
accommodated by axions via a suitable adjustment of the initial axion mis-alignment angle 

®We have also scanned a benchmark point in the hyperbohc branch/focus point region [pst of mSUGRA, 
again with a standard underabundance of neutralino dark matter. The vs. fa results look qualita- 

tively much like Fig. bl We do not present these results here in the interests of brevity. 
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10^ i(y* 10^ 10^ 10^ lo** 10" 10" 
rris (GeV) 



Figure 6: Points with i^^i^^ ^ ^-^^ obtained through the random scan described in the text. For 
comparison, we also show the curves for T£, < 5 MeV and T£, > Te and 17- ft.^ = 0.11 obtained 



from Eq. 3.8 



In Fig. 10, we show the required value of 9i needed to enforce Cl^ h? + VLah? = 0.11. 



At low fa, the points satisfying ^^Ji^ < 0.11 have axinos and saxions decaying before the 
neutralino freezes out and consequently before axions start to oscillate. Hence the axion 
relic density is not affected by the entropy injection of axinos/saxions and is given by the 
standard expression [p9[: 

/ f \ 7/6 

From the above equation, we see that as fa increases, 9i must decrease in order to maintain 
^Zi^^ + Oq/i^ = 0.11. This behavior is clearly seen in Fig. |lO| for fa < 10^^ GeV. Once 
fa becomes sufficiently large so axinos and saxions start to decay after the axion starts to 
oscillate, the entropy injected from saxions and axinos considerably dilute the axion relic 
density, thus allowing for larger 6^ values. However, as seen in Fig. [l^, this only happens 
for the BBN-forbidden solutions at fa ^ 10^^ GeV. The only BBN-allowed points at large 
fa with < 0.11 are the ones where the saxion production is either suppressed or 

where it decays before neutralino freeze-out. In this case there is no significant entropy 
injection and the axion relic density is once again given by Eq. 3.10| . Thus, extremely 



small values of 9i are required in order to suppress the axion relic density at large fa, as 
seen in Fig. |l^. Therefore these points tend to have neutralino domination of the dark 
matter density, rather than axion domination. For these points, the large neutralino halo- 
annihilation rates, enhanced by the yl-resonance, may lead to visible production rates of 
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Figure 7: Evolution of radiation, neutralino, axion, saxion, axino and gravitino energy densities 
versus scale factor R starting at T — Tr. We adopt an mSUGRA SUSY model with parameters 
(too, mi/2, Ao, tan/?, si5n(/i)) — (400 GeV,400 GeV, 0,10, +). The PQ parameters are listed on the 
plot. 



7S, e'^s and ps in cosmic ray detectors |£6[, while corresponding direct neutralino detection 
rates may remain low. 

3.5 Benchmark BM4: AMSB with wino-hke neutrahno 



In Fig. 11, we plot i}^ h for an anomaly- mediated SUSY breaking model (AMSB) with 
a wino-like neutralino ||4ll| . We choose the gaugino-AMSB model with mo ~ ~ 0, since 
this model avoids tachyonic sleptons without introduction of an additional scalar mass 
parameter [42 1 . Model parameters are (?Ti3/2, tan/3, sign{^)) = (50 TeV, 10, +), with a 
standard abundance il-'^/i^ ~ 0.0016, far below the measured value. From the figure, we 
see that for ~ 10^ — 10^^ GeV, the neutralino abundance can be enhanced and brought 
into accord with measured values. For low /„, axino production and decay augments the 
abundance, while for high /„, saxion production and decay both augments and dilutes 
the abundance. In this case, as with BM3, the standard underabundance of DM can be 
augmented and brought into accord with cosmological measurements. Unlike BM3, there 
exists no intermediate range of fa which is always excluded by the production of too much 
neutralino DM. The PQ scale can be as large as fa ~ 10^^ GeV, in accord with expectations 
from string theory. 

Originally, Moroi and Randall had proposed augmenting the relic wino abundance 
from AMSB via moduli production and decay |]57|, ^ 59|. Here, we see that an alternative 
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Figure 8: Calculated neutralino relic abundance versus fa from mSUGRA SUSY model BM2. The 
spread in dots is due to a scan over PQ parameters fa, Tr, rria, ms, 9s- 




Figure 9: Calculated neutralino relic abundance versus fa from mSUGRA SUSY model BM3. The 
spread in dots is due to a scan over PQ parameters fa, Tr, rria, ms, Og- 



mechanism introducing the several PQMSSM fields can also do the job. Direct and indirect 
detection rates for wino-like neutralinos have been presented in Ref. BQ]. 



4. The case of very large 6s and < 2mg 

From the results presented in Sees. 3.2-3^, it seems difficult to suppress the neutralino 
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Figure 10: Values of the axion mis-alignment angle 9i for the points in Fig. 
The parameter 6i is chosen such as ^^Ji^ + ^ah^ = 0.11. 



I with n^h"^ < 0.11. 



inoAMSB (m , m . A . tan^. \l) = (0 GeV. 161.4 Ge V, 0. 10. > 0) 




fa (GeV) 



Figure 11: Calculated neutralino relic abundance versus fa from inoAMSB model BM4. The 
spread in dots is due to a scan over PQ parameters fa, Tji, rria, iris, Qs- 



CDM abundance below the standard neutralino abundance. This conclusion relies on the 
fact that the (CO) saxion production and decay are correlated through the value of the PQ 
scale, since the saxion field strenght (s(x) = 6sfa)~ which sets the amplitude of the coherent 
oscillations- is assumed to be of order fa {6s ~ 0(0.1 — 10)). Hence large (CO) saxion 
production only happens at large fa values and leads to late decaying saxions, usually 
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violating the BBN bounds. As also discussed above, the BBN bounds can be avoided if 
saxions have masses in the multi-TeV range, but then the s ^ gg decay is kinematically 
allowed and the injection of neutralinos enhances the CDM abundance. However, if the 
saxion field strength (s(x)) is not set by the PQ breaking scale, but by a much larger 
scale, such as the reduced Planck mass (as suggested in some models ||5^), it is possible 
to envision a large production of coherent oscillating saxions even at small fa values. 
In this scenario, small fa easily satisfies the BBN bounds, allowing for sub-TeV saxion 
masses, such as rUg < 2mg. Thus, assuming s{x) ^> /„ {6s ^ 1), it is possible to have 
large saxion production via coherent oscillations, small fa values and small saxion masses 
without violating the BBN bounds. In this case, if nis < 2mg, saxion decay leads to large 
entropy production, but does not inject neutralinos. 

To illustrate the large dg (3> 1) scenario, in Fig. ^ we fix the initial saxion field 
strength to s{x) = Ogfa = 5 x 10^'' GeV, but allow /„ to vary and compute the neutralino 
and axion CDM abundances assuming = rus = rria = 1 TeV, Tr, = 10^ GeV, 6i = 0.5 
and the BMl benchmark point. In this case, < 2mg so that if saxions can dominate 
the energy density of the universe, they only lead to entropy dilution, and not CDM 
production. From the plot, we see that for low fa the neutralino abundance is enhanced 
due to large thermal production of axinos and their decay to neutralinos. As fa increases, 
thermal production of axinos and saxions becomes suppressed, while the saxion decay 
temperature decreases, leading to increased entropy dilution of the neutralino abundance. 
At fa ~ 10^^ GeV, ^^_^h? drops below 0.1, and the BMl point becomes allowed in the 
PQMSSM. Meanwhile, the axion abundance will also suffer entropy dilution as fa increases, 
but this is counterbalanced by an increasing axion field strength, which leads to greater 
axion production via COs: the net result is an almost flat value of O^/i^ as fa varies. Once 
fa increases past ~ 10^3 QgY^ i^i^g 

saxion becomes sufficiently long-lived that the model 
begins to violate BBN bounds. While this scenario does provide a strong dilution of dark 
matter relics, we note here that the values of Og needed are in the range 9s ~ 10^ - 10^ so 
that the saxion field strength is far beyond the value of fa and must be given by another 
physics scale. 

5. Conclusions 

In this paper we have presented the results of a calculation of mixed axion/neutralino CDM 
abundance using a set of eight coupled Boltzmann equations. The calculation improves 
upon previous results in several respects: 1. it allows for non-constant values of {crv), 
as occurs for bino-like neutralinos, where s-wave annihilation is suppressed, 2. it allows 
for interplay between neutralino enhancement via axino production and decay, while si- 
multaneously allowing for neutralino production and dilution via saxion production and 
decay, 3. it includes the effect of gravitino production and decay (not a big effect for the 
parameters presented here) and 4. it moves out of the "sudden decay" approximation and 
allows for continuous axino, saxion and gravitino decay. Our calculation allows for the 
accurate estimate of mixed axion/neutralino abundance for general choices of PQMSSM 
parameters. 
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Figure 12: Calculated neutralino and axion abundance versus fa for SUSY model BMl with 6*^/, 
fixed at 5 X 10" GeV, and with = ttIs ^ nia 1 TeV, Tr = fO^ GeV and Oi = 0.5. 



In most gravity-mediated SUSY breaking models with gaugino mass unification, it is 
typically the case that the lightest SUSY particle is a bino-like neutralino. Over most of pa- 
rameter space of models such as mSUGRA, bino-like neutralinos give rise to a dark matter 
abundance far above WMAP limits|l^, and hence vast regions of parameter space are con- 
sidered as excluded due to overproduction of neutralino dark matter. In this paper, we have 
shown that if the MSSM is extended to the PQMSSM- including an axion/saxion/axino 
supermultiplet- then SUSY models with a standard overabundance of neutralinos are typ- 
ically still excluded, even for very large values of fa ^ 10^^ — 10^^ GeV, where it might be 
expected that a high rate of entropy production from saxion decay would dilute the DM 
abundance. Here, we find that s ^ gg compensates against entropy dilution, and prevents 
the neutralino abundance from dropping into the measured range, unless the saxion decays 
are in violation of BBN bounds on late-decaying neutral particles. As noted earlier, our 
conclusion depends on at least three assumptions. First, we implemented the standard 
thermal axino production rates as calculated in the Ref'spl], p9| , These rates should 
apply in supersymmetric versions of the KSVZ model where PQ-charged matter multiplets 
$ exist at or around the PQ breaking scale fa- In a recent publication pI| , it has been 
shown that in the SUSY DFSZ model, thermal axino production rates can be enhanced or 
diminished compared to their KSVZ values depending on PQMSSM parameters. Secondly, 
we assumed that saxion decay is dominated by two-body modes into gluon and gluino pairs. 
In the DFSZ model, decays into Higgs pairs or aa may also contribute, and even dominate 
the saxion decay modes. Thirdly, we have assumed saxion field strength s{x) = 6sfa is of 
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order the PQ-breaking scale fa, i.e. that 0^ ~ 1. We have also shown in Sec. ^ that if 
9s ^ 1 and < 2mg, then CO-produced saxions can dominate the universe and dilute 
all thermal relics while avoiding BBN constraints. 

In the case of a standard underabundance of neutralino CDM, a wide range of fa 
values are permitted, and can augment the neutralino DM into the measured range. In 
cases where the neutralinos still maintain an underabundance, the remaining abundance 
can be accommodated by axions. In these cases of a standard underabundance of neutralino 
DM, the PQ scale fa can be pushed into the 10^"^ — 10^^ GeV range, which is closer to 
expectations from string theory. For the case of very high fa, then we typically expect the 
DM to be neutralino rather than axion dominated, since the neutralino abundance cannot 
be suppressed too much without violating BBN constraints. 
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A. Boltzmann Equations for the PQMSSM 

As discussed in Sec. H we assume the following set of coupled differential equations: 



Hi 



-SHm - Tirrii^ + [(nf (T))^ - nf]M, + BR{j, i)r,m,-^ , 

Pi j Pj 



S =^Y^ BR{i, X)rimini , (A.l) 

i 



with H given by: 



where pr is the total energy density. 

In order to simplify the above equations we define: 

X = ln{R/Ro), Ni = ln(ni/so), and Ns = HS/So) (A.3) 
so we can write Eq's. [A.l| as: 

Ns = exp[iV, + 3x- Ns] (A.4) 

i 

3 _ T m ^ ^g^, , ^ Mi„.|(<!)^ i| (A,5) 

B Pi/rii H pj/uj Ui H \ni J 

where ' = d/dx and rii is given by Ui = sqc^'. 

The above equation for Ni also applies for coherent oscillating fields, if we define: 

Ni = ln(nj/so), and rn = pi/rrii (A.6) 
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so 



N' 



(A.7) 



where we assume that the coherent oscillating component does not couple to any of the 
other fields. 

Since H depends on the energy densities, to solve the above equations we must compute 
Pi from rij. However, even for particles following a thermal distribution, the energy density 
for each component cannot be directly obtained from n,, unless the chemical potential 
(Hi) is also given. Nonetheless, /Xj(T) is usually small in the relativistic regime, while in 
the non-relativistic regime we always have pi = rriirii. Therefore, assuming that the fields 
follow a thermal distribution, a good approximation for pi as a function of Ui is given by: 



ra,- X < 



m 

^^K2{mjT,) + '^^^ 
^^{3) 30 



if Ti < rrii/lO 

if mi/10 <Ti< 3mi/2 

if 3mi/2 < Ti 



where the modified Bessel functions, Ki and K2, are necessary to describe a smooth 
relativistic/non-relativistic transition and Np = 1(7/6) for bosons (fermions). 

The only remaining piece of information necessary for computing pi and H and solv- 
ing the Boltzmann equations is the definition of temperature for each component. The 
radiation temperature can be directly obtained from Ns and x: 

T=(^^^y\neMNs/3-x]. (A.9) 

For thermal fluids in equilibrium we always have Tj = T, but once they decouple, this is 
no longer true. However, the temperature of relativistic fluids scales as T oc R~^, while 
non-relativistic fluids have T oc R^^. Thus, we approximate Tj by 



X < 



T 



dec 



, if coupled 

, if Tj > ?)rail2 and decoupled 



|mj , if Tj < 2>mi/2 and decoupled 



(A.IO) 



where T^'"', Rf" and Rf^ are the decoupling (freeze-out) temperature, the scale factor at 
freeze-out and the scale factor at the non-relativistic transition (Tj = 3m-j/2), respectively. 
If the fluid was never in thermal equilibrium, we take Tf'^'^ = Tji. For coherent oscillating 
fluids we always have Tj = 0^. 

Pi (Eq. 



Eq's. [A.4| and [A.5| , with the auxihary equations for H (Eq. [A.2| ), pi (Eq. [A.8| ) and Tj 
(Eq. |A.10 ) form a set of closed equations, which can be solved once the initial conditions 
for the number densities (nj) and entropy {S) are given. The initial entropy Sq is trivially 
obtained, once we assume a radiation dominated universe at T = Tr: 



S{Tr) = —g4Tn)nRl 



(A.ll) 



^In principle, the approximations in Eq's. A. 8 and A.IO can be avoided if we include equations for the 



chemical potentials fii{T). However, for simplicity, we use Eq's. [A.8| and A.IC instead 
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For thermal fluids we take the initial number density as 



, if {av)in':yH\T=Tn < 10 

n1%Tn) , if {av)in:iyH\T=Tn > 10 



(A.12) 



while for coherent oscillating fluids the initial condition is set at the beginning of oscilla- 
tions: 



n^{Tt 



(A.13) 



where T°^'^ is the oscillation temperature, given by 3H{T°^'^) = mi{T°'^^) and the initial 
energy density for oscillations. For the oscillating saxion and axion [^9| fields the initial 
energy densities are given by: 



1.44- 



mm 



2.1 X 10"^ 



27r2g,(Tfi)r| 
45 



10^ 



10 



12 



where f{9i) = ln[e/(l — ^^/vr^)] and 9ifa and Ogfa are the initial axion and saxion field 
amplitudes. The definition of accounts for the possibility of saxion oscillations beginning 
during inflation (if Tr < Tosc)- 

In order to compute the source term in Eq. |A.1| , we must specify the annihilation 
cross-sections {crv)i, the branching ratios BR{i,j) and BR{i,X) and the the decay widths 
Fj. The annihilation cross-sections for axions, saxions, axinos and gravitinos are given by 
the expressions |29, 62 1 

n6 



{<yv)a 



10" 
10" 



^ [4.191n(1.5/52) + l.68x0(r-ma(T))] 
3 3.87 X 9{T-ma 



5 at 



(fa) 



Mg = — ^ X 



Ml 



72gl\n{l.27l/gs){l + 




3m1' 
G 



23.863c/70 '^ - 0.784 , if c/s > 0.35 

4.47 + 31 ln(1.4/5^) - 0.784 , if 5^ < 0.35 



Ml 



+ llg'^ ln(1.266/5')(l 



3m?- ^ 

G 



while {av)^^ (T) is extracted from IsaReD [^] . The second term in the expressions for {av)a 
and {uv)a represent contributions from 1 — )• 2 decays of particles with thermal masses. 
Therefore, these terms should not be included unless T > ma,a, as indicated by the 6 
functions above. The expression for the axino effective cross-section is set to reproduce the 
numerical results in |3C]. Since the saxion thermal production has not been computed, we 
approximate it by the axion expression: 



{av)s = {(Jv)a 



(A.14) 



with nin 
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For obtaining the various unstable particle widths, we calculate Ta from the a — >• gg, 
Zij and ZiZ partial widths as presented in Ref. [26|. For gravitino decays, we adopt the 



gravitino widths as presented in Ref. |33|. For the saxion width, we include Tg from the 
s ^ gg and s ^ gg decays as presented in Ref. ||5^. We note here that in the DFSZ 
model, it is also possible to have s ^ hh decays and possibly s — )■ aa decays. We neglect 
these latter two cases, so that our results apply to the supersymmetrized KSVZ model, 
where the gg and gg final states should dominate. 

Once the total and partial widths are known, we can easily compute the required 
branching ratios: 



BR{a, Zi) = 1, BR{s, Zi) = 2 x 



r(s ^ gg) 



BR{G,Zi) = 1 



(A.15) 



The factor 2 in BR{s,Zi) takes care of the multiplicity of neutralinos for each saxion 
cascade decay. While the s ^ gg decay width is always dominant, we showed in Sec. ^ 
that s ^ gg plays a crucial role in the PQMSSM dark matter cosmology. 

Finally, we assume that the branching ratios for computing the energy injection into 
the thermal bath from unstable particle decays are given by: 



BR{a,X) = BR{s,X) = BR{G,X) = 1. 



(A.16) 



Although some of the decay energy is lost into neutralinos (except for s ^ gg decays), 
we assume that in the final product of the cascade decay of axinos, saxions and gravitinos 



most of the initial energy has been converted into radiation, so Eq. A.16 consists in a good 
approximation. 
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